library(raster)
library(rgdal)
library(rgeos)
library(plyr)
library(reshape)
library(reshape2)
library(maps)
library(maptools)
library(foreign)


## Railroads 
## Data source: Jeremy Atack, Historical Geographic Information Systems (GIS) database of U.S. Railroads
## https://my.vanderbilt.edu/jeremyatack/data-downloads/
## Go to "Zipped SHP files" > Download "Railroads, 1826-1911"
rail <- readOGR("RR1826-1911Modified0509161", "RR1826-1911Modified050916")
rail1850 <- rail[rail$InOpBy < 1850,]
rail1860 <- rail[rail$InOpBy < 1860,]

## State boundaries data can be obtained from here:
## https://digital.newberry.org/ahcb/downloads/united_states.html
## Go to "US Historical States & Territories (1783-2000)" > Generalized .001 deg. tolerance	> GIS Files (14MB)
histstates <- readOGR("US_AtlasHCB_StateTerr_Gen001/US_HistStateTerr_Gen001_Shapefile", "US_HistStateTerr_Gen001")
histstates <- spTransform(histstates, CRS(proj4string(rail)))
states1860 <- histstates[histstates$START_N <= 18600101 & histstates$END_N >= 18600101,]
states1860 <- states1860[states1860$ABBR_NAME %in% c("FL", "GA", "LA", "MD", "NC", "SC", "VA", "AL", "AR", "KY", "MS", "MO", "TN", "TX"),]

# Overlay maps
intersect1850 <- raster::intersect(rail1850, states1860)
intersect1860 <- raster::intersect(rail1860, states1860)

## County boundaries can be obtained from here:
## https://digital.newberry.org/ahcb/downloads/united_states.html
## Go to: US Historical Counties (1629-2000) > Generalized .001 deg. tolerance > GIS Files (52MB)
histcounties <- readOGR("US_AtlasHCB_Counties_Gen001/US_HistCounties_Gen001_Shapefile", layer = "US_HistCounties_Gen001")
histcounties$fips <- as.numeric(as.character(histcounties$FIPS))
histcounties$fips[which(histcounties$NAME == "CARROLL (ext)")] <- 22035
histcounties <- spTransform(histcounties, CRS(proj4string(rail)))

# 1860 counties
data("county.fips")
cc1860 <- which(histcounties$START_N <= 18590101 & histcounties$END_N >= 18590101)
counties1860 <- histcounties[cc1860,]
kk <- as.character(counties1860$NAME[counties1860$STATE_TERR == "South Carolina"])
kk <- sub(" Dist.", "", kk)
kk <- paste("south carolina,",tolower(kk), sep = "")
kk.fips <- NA
kk.fips[kk %in% county.fips$polyname] <- unlist(lapply(kk, function(x) county.fips$fips[county.fips$polyname == x]))
counties1860$fips[counties1860$STATE_TERR == "South Carolina"] <- kk.fips
counties1860 <- counties1860[!is.na(counties1860$fips),]
newpolys <- unionSpatialPolygons(counties1860, counties1860$fips)
kk <- unique(as.data.frame(counties1860)[,c("fips", "STATE_TERR")])
rownames(kk) <- as.character(kk$fips)
counties1860 <- SpatialPolygonsDataFrame(newpolys, data = kk)

# 1860 Census data
census1860 <- read.csv("census1860.csv")
counties1860 <- merge(counties1860, census1860, by = "fips")

### Figure A1: Railroads and County Slave Share (% Pop), 1860

pslave.bins <- cut(counties1860$sprop, breaks = seq(0,1, length = 6))
mcols <- c("#eff3ff", "#c6dbef", "#bdd7e7", "#6baed6", "#3182bd", "#08519c")
leg.text <- c("0-20%", "21-40%", "41-60%", "61-80%", "81-100%")
names(mcols) <- levels(pslave.bins)

southern.states <- c("Alabama", "Arkansas", "Florida", "Georgia", "Kentucky", "Louisiana",
                     "Mississippi", "Missouri", "North Carolina", "South Carolina", "Tennessee", "Texas",
                     "Virginia", "West Virginia", "Maryland")

# counties
plot(counties1860[counties1860$STATE_TERR %in% southern.states,], col = mcols[pslave.bins[counties1860$STATE_TERR %in% southern.states]], border = NA)
plot(states1860[states1860$FULL_NAME %in% southern.states,], border = "grey30", add = TRUE, lwd = 0.25)
legend("topleft", horiz = FALSE, col = mcols, legend =  leg.text, bty = "n", fill = mcols, border = "white")
# railroads
plot(intersect1860, add=TRUE, col='red', lwd=1)
plot(intersect1850, add=TRUE, col='black', lwd=1)
legend("bottomleft", horiz = FALSE, col = c("black", "red"), lty=1:1, legend = c("1850", "1860"), 
       bty = "n", title = "Railroads", border = "white")



### Figure A2: Apportionment Status and Enslaved Share (% of Pop.), 1860

plot(states1860[states1860$FULL_NAME %in% c("Florida", "Georgia", "Louisiana", "South Carolina",
                                            "North Carolina", "Virginia", "West Virginia", "Maryland",
                                            "Kentucky", "Missouri", "Arkansas", "Tennessee", "Texas",
                                            "Alabama", "Mississippi"),], border = "grey10", lwd = 0.25)

plot(states1860[states1860$FULL_NAME %in% c("Florida", "Georgia", "Louisiana", "South Carolina"),], col = "#fb6a4a", add = TRUE, border = NA)
plot(states1860[states1860$FULL_NAME %in% c("North Carolina", "Virginia", "West Virginia"),], col = "#fcbba1", border = NA, add = TRUE)
plot(states1860[states1860$FULL_NAME %in% c("Maryland"),], col = "#fee5d9", border = NA, add = TRUE)

plot(states1860[states1860$FULL_NAME %in% c("Kentucky", "Missouri"),], col = "#deebf7", border = NA, add = TRUE)
plot(states1860[states1860$FULL_NAME %in% c("Arkansas", "Tennessee", "Texas"),], col = "#9ecae1", border = NA, add = TRUE)
plot(states1860[states1860$FULL_NAME %in% c("Alabama", "Mississippi"),], col = "#3182bd", border = NA, add = TRUE)

mcols1 <- c("#fb6a4a", "#fcbba1", "#fee5d9")
mcols2 <- c("#3182bd", "#9ecae1", "#deebf7")
leg.text1 <- c("> 40%", "25-33%", "0-20%") 
leg.text2 <- c("> 40%", "25-33%", "0-20%")
legend("topleft", horiz = FALSE, col = c("#fb6a4a", "#fcbba1", "#fee5d9"), title="Slave share\nMalapportioned States", y.intersp=.6, legend= leg.text1, bty = "n", fill = mcols1, border = "white")
legend("bottomleft", inset= -.05, horiz = FALSE, col = c("#3182bd", "#9ecae1", "#deebf7"), title="Non-malapportioned States", y.intersp=.6, legend= leg.text2, bty = "n", fill = mcols2, border = "white")



